import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
import pandas as pd

def sigmoid(x):
    return 1/(1+np.exp(-x))
def dsigmoid(x):
    return x*(1-x)                                                      #定义sigmoid的导数函数
def update(v,w,b1,a,b2,n,v1,v2,t):                                      #定义迭代函数，v的位置对应Iw,w的位置对应Hw下面分别用第一个和第二个矩阵表示,n表示迭代数据笔数，t是目标值
    global O
    for i in range(n):
        xw=np.array([x[i]])                                             #增加第i笔数据的维度
        tw=np.array([t[i]])                                             #增加第i笔目标值的维度
        h=sigmoid(np.dot(xw,v))                                         #第i笔数据乘第一个矩阵然后sigmoid,的到隐藏元的输出
        H=np.array([np.insert(h,0,1)])                                  #添加h0=1，因为insert让维度下降，再升维                             
        O=np.dot(H,w)                                                   #隐藏元的输入乘第二个矩阵得到输出元的输出
        H_delta=(tw-O)                                                  #计算第二矩阵的delta,是目标矩阵减输出矩阵
        w1=np.delete(w,0,0)                                             #删除第二个矩阵的第一列(h0不需要返回）
        I_delta=np.dot(H_delta,w1.T)*dsigmoid(h)                        #计算第一个矩阵的delta,其值等于第二个矩阵的delta与w1.T矩阵相乘，再点乘隐藏元输出的dsigmoid
        wc=b2*np.dot(H.T,H_delta)                                       #第二个矩阵的权值调整量为对应学习率乘隐藏元输出的转置与第二个矩阵的delta相乘
        vc=b1*np.dot(xw.T,I_delta)                                      #第一个矩阵的权值调整量为对应学习率乘i笔数据的转置与第一个矩阵的delta相乘
        w=w+a*v2+wc                                                     #第二个矩阵的调整为上一次矩阵加惯性率乘上一次矩阵改变量加权值调整量
        v=v+a*v1+vc                                                     #第一个矩阵的调整为上一次矩阵加惯性率乘上一次矩阵改变量加权值调整量
        v2=wc+a*v2                                                      #第二个矩阵改变量赋值给v2
        v1=vc+a*v1                                                      #第一个矩阵的改变量赋值给v1
    return w,v,v1,v2,O
def cycle(a,b1,b2,t):
    global Hw,Iw,v1,v2,oL,O1
    fig=plt.figure()
    ax1=fig.add_subplot(141)                                            
    ax1.set_xlabel('x1')
    plt.axis([0,1,-1,1])
    ax2=fig.add_subplot(142)
    ax2.set_xlabel('x1')
    plt.axis([0,1,-0.5,2.5])
    ax3=fig.add_subplot(143)
    ax3.set_xlabel('x1')
    plt.axis([0,1,-2,2])
    ax4=fig.add_subplot(144)
    ax4.set_xlabel('x1')
    plt.axis([0,1,-1.5,3.5])     
    ax1.plot(x1[:,1],t[:,0],c='black')
    ax2.plot(x1[:,1],t[:,1],c='black')
    ax3.plot(x1[:,1],t[:,2],c='black')
    ax4.plot(x1[:,1],t[:,3],c='black')
    plt.ion()

    In=1
    Iw=np.random.random((In+1,H1))                                          #依据输入和隐藏元随机生成初始矩阵，矩阵的行是输入元加1，列是隐藏元
    Hw=np.random.random((H1+1,4))                                         ##依据隐藏元和输出元随机生成初始矩阵，矩阵的行是隐藏元加1，列是输出元
    v1=np.zeros((In+1,H1))                                                  #生成对应的0矩阵，作为第一次权值迭代的上一次改变量
    v2=np.zeros((H1+1,4))
    for i in range(Iter):
        Hw,Iw,v1,v2,O1=update(Iw,Hw,a,b1,b2,xm.shape[0],v1,v2,t)       #每次循环输出的Hw,Iw,v1,v2为下一次迭代的输入，update的数据笔数为输入数据的笔数
        h_p=sigmoid(np.dot(x,Iw))                                           #进行一次循环计算一次1001笔data的隐藏层输出
        h0_p=np.ones((xm.shape[0],1))                                       #添加h0          
        H_p=np.hstack((h0_p,h_p))                                           #计算一次1001笔data的输出output
        OL=np.dot( H_p,Hw)                      
        if(i/(Iter-1)==1):
            print('Final Weight_In_Hid:\n',Iw)                              #打印循环的最后一次的权值矩阵
            print('Final Weight_Hid_Out:\n',Hw)
        ax1.set_title('Iteration=%s'%(i+1)
        if(I==1):                                                           #是否要动态演示的信号   
            plt.pause(0.1)
            try:                          
                ax1.lines.remove(line1[0])
                ax2.lines.remove(line2[0])
                ax3.lines.remove(line3[0])
                ax4.lines.remove(line4[0])
            except Exception:
                pass                                                        
            if(In==1):
                line1=ax1.plot(x[:,1],OL[:,0],c='red')
                line2=ax2.plot(x[:,1],OL[:,1],c='red')
                line3=ax3.plot(x[:,1],OL[:,2],c='red')
                line4=ax4.plot(x[:,1],OL[:,3],c='red')
    if(I==0):                                                    #不动态演示就画x1,最后一个output对应的8张图
        ax1.plot(x1[:,1],OL[:,0],c='red')
        ax2.plot(x1[:,1],OL[:,1],c='red')
        ax3.plot(x1[:,1],OL[:,2],c='red')
        ax4.plot(x1[:,1],OL[:,3],c='red')  
    plt.ioff()
    plt.show()
      
if __name__=='__main__':
    path=('./data.xlsx')
    Data1=pd.read_excel(path,0,index_col=0)                         #输入元和输出元为1时打开第一张表
    xm=np.array(Data1)                                                         
    x0=np.ones((xm.shape[0],1))
    x1=np.hstack((x0,xm))                                               #水平合并对应x0
    x=x1[:,:2]                                                          #进行迭代数据是前两列                                                      
    t=x1[:,2:]
    print('='*30,'\tBackpropogaration Network\t','='*15)
    print('Artificial Network with Backpropogaration learning\n','='*80)
    
    num=int(input('Best Parameters <0-No, 1-Yes>：'))                        #是否要用设计好的的参数
    if(num==0):                                                              #否的话手动输入各项参数
        H1=int(input('Hidden layer neurons：'))                              #隐藏元个数赋值给H1
        Iter=int(input('Total Iteration:'))                                  #迭代次数赋值给Iter
        a=float(input('Inertia rate<range in [0,1]>:'))                      #惯性率赋值给a，两个惯性率都为a
        b1=float(input('Learning rate in Weight_In_Hid<range in [0,1]>:'))   #Weight_In_Hid的学习率赋值给b1
        b2=float(input('Learning rate in Weight_Hid_Out<range in [0,1]>:'))  #Weight_Hid_Out的学习率赋值给b2
    if(num==1):                                                              #是的话用多次实验得到的参数
        H1=20
        Iter=999
        a,b1,b2=0.1,0.1,0.001
    I=int(input('Drawing per iteration <0-No,1-Yes>:'))         #输入1就实现第一张到最后一张的动画演示，In=0时没有动画
    cycle(a,b1,b2,t)           #迭代
